estimate_matches<-function(block_vars, data){
    df_weights = data[,which(names(data) %in% block_vars)]
    df_weights = data.frame(C=data$C, df_weights)[data$T==1,]

    names(df_weights)[-1]<-block_vars
    model = glm(C~., data=df_weights, family='binomial')
    weights = predict(model, data, 'response')
    psw_full = mean(data$Y[data$T == 1 & data$C == 1]) - sum(data$Y[data$T==0]*weights[data$T==0])/sum(weights[data$T==0])
    varY1 = var(data$Y[data$T == 1 & data$C == 1])/sum(data$C)
    varY0 = var(data$Y[data$T==0])*sum(weights[data$T==0]^2)/(sum(weights[data$T==0])^2)
    psw_full_se = sqrt(varY1+varY0)

    #Set up matched data frame
    matching=Match(Y=data$Y, Tr=data$T, X=data[,which(names(data) %in% block_vars)], replace=FALSE, M=1)
    match_index = 1:length(unique(matching$MatchLoopC[,1]))

    data$match_index = NA
    data$match_index[matching$MatchLoopC[,1]]<-match_index
    data$match_index[matching$MatchLoopC[,2]]<-match_index
    df_match<-na.omit(data)

    complier_matches<-df_match$match_index[which(df_match$C==1)]

    #BLOCK IV:
    model_iv_block = iv_robust(Y~C|T, data=df_match, fixed_effects=~match_index,
                           try_cholesky=TRUE, se_type="HC1")
    iv_block = coef(model_iv_block)[1]
    iv_block_se = model_iv_block$std.error[1]

    #BLOCK DIM:
    df_block = df_match %>%
        group_by(match_index) %>%
            summarize(
                block_weight = sum(C[T==1])/sum(df_match$C),
                mdim=mean(na.omit(Y[T==1]))-mean(na.omit(Y[T==0])),
                mY = mean(na.omit(Y[T==0]))
            )

    bdim = sum(df_block$block_weight*df_block$mdim)
    bdim_se = sqrt(1/(sum(df_match$C)-1) * sum(df_block$block_weight*(df_block$mdim - mean(df_block$mdim))^2))

    #BLOCK PSW
    df_weights = df_match[,which(names(df_match) %in% block_vars)]
    df_weights = data.frame(C=df_match$C, df_weights)[df_match$T==1,]
    names(df_weights)[-1]<-block_vars
    model = glm(C~., data=df_weights, family='binomial')
    weights = predict(model, df_match, 'response')
    weights_block = ifelse(df_match$T == 0, weights, ifelse(df_match$C==1, 1, 0))
    model_psw_block = lm_robust(Y~T, weights=weights_block, clusters=match_index,
        data=df_match)
    psw_block = coef(model_psw_block)[2]
    psw_block_se = model_psw_block$std.error[2]

    #SET UP SRS:
    df_treatment = data[data$T==1,]
    df_control = data[data$T==0,]
    n_treatment = sum(df_match$T)
    n_control = sum(df_match$T==0)
    df_srs = rbind(
        df_treatment[sample(1:nrow(df_treatment), n_treatment, replace=FALSE),],
        df_control[sample(1:nrow(df_control), n_control, replace=FALSE),])

    #PSW (SRS)
    df_weights = df_srs[,which(names(df_srs) %in% block_vars)]
    df_weights = data.frame(C=df_srs$C, df_weights)[df_srs$T==1,]

    names(df_weights)[-1]<-block_vars
    model = glm(C~., data=df_weights, family='binomial')
    weights = predict(model, df_srs, 'response')
    psw = mean(df_srs$Y[df_srs$T == 1 & df_srs$C == 1]) - sum(df_srs$Y[df_srs$T==0]*weights[df_srs$T==0])/sum(weights[df_srs$T==0])
    varY1 = var(df_srs$Y[df_srs$T == 1 & df_srs$C == 1])/sum(df_srs$C)
    varY0 = var(df_srs$Y[df_srs$T==0])*sum(weights[df_srs$T==0]^2)/(sum(weights[df_srs$T==0])^2)
    psw_se = sqrt(varY1+varY0)


    #IV (SRS)
    model_iv = iv_robust(Y ~ C | T , data=df_srs, try_cholesky=TRUE, se_type="HC1")
    iv = coef(model_iv)[2]
    iv_se = model_iv$std.error[2]

    at = mean(df_srs$Y[df_srs$C==1]) - mean(df_srs$Y[df_srs$T==0])

    return(list(
        summary = data.frame(at, iv, iv_se, psw, psw_se, iv_block, iv_block_se,
            bdim, bdim_se, psw_block, psw_block_se,
            psw_full, psw_full_se),
        matching_results = matching,
        n_units = nrow(df_match),
        revealed_compl = data.frame(match = sum(df_match$C)/sum(df_match$T), srs = sum(df_srs$C)/sum(df_srs$T))
    ))
}